Recommendation for reading and cutoffs: https://www.pnas.org/doi/10.1073/pnas.2104429118
Overview: 1. PCoA Analysis - complete microbiome using Unifrac? distances 2. Alpha diversity and betadispersion 3. Determine the core across control samples 4. How much of the community is captured? 5. How does community change with exposure to ABS? 6. Which core members are differential abundant at TP4,TP5,TP6?
library(tidyr)
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.3 ✔ purrr 1.0.2
✔ forcats 1.0.0 ✔ readr 2.1.4
✔ ggplot2 3.4.4 ✔ stringr 1.5.0
✔ lubridate 1.9.3 ✔ tibble 3.2.1── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(ggplot2)
library(vegan)
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-4
library(phyloseq)
library(ggrepel)
library(patchwork)
library(microbiome)
microbiome R package (microbiome.github.com)
Copyright (C) 2011-2022 Leo Lahti,
Sudarshan Shetty et al. <microbiome.github.io>
Attaching package: ‘microbiome’
The following object is masked from ‘package:vegan’:
diversity
The following object is masked from ‘package:ggplot2’:
alpha
The following object is masked from ‘package:base’:
transform
library(RColorBrewer)
plotting conditions
trt16<-c('ABS1','ABS2','Control','Experimental', 'FASW','Priming','blank', 'FAASW' )
l2p16<-c('indianred1','deeppink4','gray16','burlywood','burlywood','burlywood', 'gray81', 'darkgreen')
names(l2p16) = trt16
shpsc<-c(0,1,3,4,8,15,16)
names(shpsc)=c('T0','T1','T2','T3','T4','T5','T6')
theme_new<- theme_bw() + theme(panel.grid.major = element_line(color='gray90'), panel.grid.minor = element_line(color='gray90'), text= element_text(size = 15), axis.text.x = element_text(angle = 65,hjust = 1), legend.text = element_text(size = 20), axis.title = element_text(size = 15), title = element_text(size = 20))
read in cleaned data
ps<-readRDS(file = 'ps.cleantree.noblanks.rds')
ps<-prune_samples(!(sample_names(ps) %in% c('EB1_S94', 'EB1_S95','EB2_S95')),ps) #remove EB samples
ps
phyloseq-class experiment-level object
otu_table() OTU Table: [ 414 taxa and 82 samples ]
sample_data() Sample Data: [ 82 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 414 taxa by 7 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 414 tips and 412 internal nodes ]
refseq() DNAStringSet: [ 414 reference sequences ]
ps.red<-subset_samples(ps, Time %in% c('T0','T1','T4','T5','T6'))
p<-plot_richness(ps.red, measures=c("Observed"))
p_df<-p$data
a<-ggplot(p_df, aes(y=value, x=Condition, fill=Condition)) +
geom_boxplot(alpha=0.5, outlier.shape = NA, color='black') +
scale_fill_manual(values=l2p16)+scale_color_manual(values=l2p16)+
geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+
theme_new + ylab('Nr. ASVs') + xlab('') + guides(fill=F, group=F,color=F) + ylim(c(0,140)) +
theme(axis.text.x = element_blank(),panel.grid = element_blank(),panel.grid.minor = element_blank(),rect = element_blank())
a
stats
p_df.trt<-p_df
hist(p_df.trt$value)
p_df.trt$Time<-as.factor(p_df.trt$Time)
p_df.trt$Condition<-as.factor(p_df.trt$Condition)
aov.trt<-aov(value ~ Condition*Time, p_df.trt)
shapiro.test(residuals(aov.trt))
Shapiro-Wilk normality test
data: residuals(aov.trt)
W = 0.99104, p-value = 0.9574
summary(aov.trt)
Df Sum Sq Mean Sq F value Pr(>F)
Condition 3 7147 2382.2 8.494 0.000159 ***
Time 4 6159 1539.6 5.490 0.001199 **
Condition:Time 4 9508 2377.0 8.475 4.16e-05 ***
Residuals 42 11780 280.5
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
aov.trt
Call:
aov(formula = value ~ Condition * Time, data = p_df.trt)
Terms:
Condition Time Condition:Time Residuals
Sum of Squares 7146.709 6158.544 9507.906 11779.600
Deg. of Freedom 3 4 4 42
Residual standard error: 16.74714
8 out of 20 effects not estimable
Estimated effects may be unbalanced
par(mfrow=(c(2,2)))
plot(aov.trt)
data_summary <- group_by(p_df.trt, Condition, Time) %>%
summarise(mean=mean(value), sd=sd(value), se=sd(value)/sqrt(length(value))) %>%
arrange(desc(mean))
`summarise()` has grouped output by 'Condition'. You can override using the `.groups` argument.
print(data_summary)
tukey <- TukeyHSD(aov.trt)
print(tukey)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = value ~ Condition * Time, data = p_df.trt)
$Condition
diff lwr upr p adj
ABS2-ABS1 8.794872 -8.180495 25.77024 0.5150285
Control-ABS1 27.747253 11.937845 43.55666 0.0001624
Experimental-ABS1 21.861538 -1.712666 45.43574 0.0777548
Control-ABS2 18.952381 3.807933 34.09683 0.0090201
Experimental-ABS2 13.066667 -10.066859 36.20019 0.4401964
Experimental-Control -5.885714 -28.177738 16.40631 0.8940480
$Time
diff lwr upr p adj
T1-T0 -5.841270 -37.65865 25.9761110 0.9845071
T4-T0 -3.384334 -33.95348 27.1848171 0.9977642
T5-T0 -9.107274 -39.47103 21.2564829 0.9115463
T6-T0 -27.233455 -57.41807 2.9511624 0.0943868
T4-T1 2.456936 -18.23848 23.1523491 0.9970659
T5-T1 -3.266004 -23.65681 17.1247993 0.9907266
T6-T1 -21.392186 -41.51526 -1.2691070 0.0322421
T5-T4 -5.722940 -24.10533 12.6594456 0.9000293
T6-T4 -23.849122 -41.93408 -5.7641683 0.0045120
T6-T5 -18.126182 -35.86175 -0.3906164 0.0430011
$`Condition:Time`
diff lwr upr p adj
ABS2:T0-ABS1:T0 NA NA NA NA
Control:T0-ABS1:T0 NA NA NA NA
Experimental:T0-ABS1:T0 NA NA NA NA
ABS1:T1-ABS1:T0 NA NA NA NA
ABS2:T1-ABS1:T0 NA NA NA NA
Control:T1-ABS1:T0 NA NA NA NA
Experimental:T1-ABS1:T0 NA NA NA NA
ABS1:T4-ABS1:T0 NA NA NA NA
ABS2:T4-ABS1:T0 NA NA NA NA
Control:T4-ABS1:T0 NA NA NA NA
Experimental:T4-ABS1:T0 NA NA NA NA
ABS1:T5-ABS1:T0 NA NA NA NA
ABS2:T5-ABS1:T0 NA NA NA NA
Control:T5-ABS1:T0 NA NA NA NA
Experimental:T5-ABS1:T0 NA NA NA NA
ABS1:T6-ABS1:T0 NA NA NA NA
ABS2:T6-ABS1:T0 NA NA NA NA
Control:T6-ABS1:T0 NA NA NA NA
Experimental:T6-ABS1:T0 NA NA NA NA
Control:T0-ABS2:T0 NA NA NA NA
Experimental:T0-ABS2:T0 NA NA NA NA
ABS1:T1-ABS2:T0 NA NA NA NA
ABS2:T1-ABS2:T0 NA NA NA NA
Control:T1-ABS2:T0 NA NA NA NA
Experimental:T1-ABS2:T0 NA NA NA NA
ABS1:T4-ABS2:T0 NA NA NA NA
ABS2:T4-ABS2:T0 NA NA NA NA
Control:T4-ABS2:T0 NA NA NA NA
Experimental:T4-ABS2:T0 NA NA NA NA
ABS1:T5-ABS2:T0 NA NA NA NA
ABS2:T5-ABS2:T0 NA NA NA NA
Control:T5-ABS2:T0 NA NA NA NA
Experimental:T5-ABS2:T0 NA NA NA NA
ABS1:T6-ABS2:T0 NA NA NA NA
ABS2:T6-ABS2:T0 NA NA NA NA
Control:T6-ABS2:T0 NA NA NA NA
Experimental:T6-ABS2:T0 NA NA NA NA
Experimental:T0-Control:T0 NA NA NA NA
ABS1:T1-Control:T0 NA NA NA NA
ABS2:T1-Control:T0 NA NA NA NA
Control:T1-Control:T0 1.5 -46.805061 49.805061 1.0000000
Experimental:T1-Control:T0 -17.6 -63.788473 28.588473 0.9946379
ABS1:T4-Control:T0 1.0 -50.640283 52.640283 1.0000000
ABS2:T4-Control:T0 -16.8 -62.988473 29.388473 0.9968728
Control:T4-Control:T0 -28.2 -74.388473 17.988473 0.7202862
tkl NA NA NA NA
ABS1:T5-Control:T0 -51.4 -97.588473 -5.211527 0.0160516
ABS2:T5-Control:T0 -14.8 -60.988473 31.388473 0.9993584
Control:T5-Control:T0 -7.5 -55.805061 40.805061 1.0000000
Experimental:T5-Control:T0 NA NA NA NA
ABS1:T6-Control:T0 -51.8 -97.988473 -5.611527 0.0146500
ABS2:T6-Control:T0 -60.4 -106.588473 -14.211527 0.0018482
Control:T6-Control:T0 -16.2 -62.388473 29.988473 0.9979815
Experimental:T6-Control:T0 NA NA NA NA
ABS1:T1-Experimental:T0 NA NA NA NA
ABS2:T1-Experimental:T0 NA NA NA NA
Control:T1-Experimental:T0 NA NA NA NA
Experimental:T1-Experimental:T0 NA NA NA NA
ABS1:T4-Experimental:T0 NA NA NA NA
ABS2:T4-Experimental:T0 NA NA NA NA
Control:T4-Experimental:T0 NA NA NA NA
Experimental:T4-Experimental:T0 NA NA NA NA
ABS1:T5-Experimental:T0 NA NA NA NA
ABS2:T5-Experimental:T0 NA NA NA NA
Control:T5-Experimental:T0 NA NA NA NA
Experimental:T5-Experimental:T0 NA NA NA NA
ABS1:T6-Experimental:T0 NA NA NA NA
ABS2:T6-Experimental:T0 NA NA NA NA
Control:T6-Experimental:T0 NA NA NA NA
Experimental:T6-Experimental:T0 NA NA NA NA
ABS2:T1-ABS1:T1 NA NA NA NA
Control:T1-ABS1:T1 NA NA NA NA
Experimental:T1-ABS1:T1 NA NA NA NA
ABS1:T4-ABS1:T1 NA NA NA NA
ABS2:T4-ABS1:T1 NA NA NA NA
Control:T4-ABS1:T1 NA NA NA NA
Experimental:T4-ABS1:T1 NA NA NA NA
ABS1:T5-ABS1:T1 NA NA NA NA
ABS2:T5-ABS1:T1 NA NA NA NA
Control:T5-ABS1:T1 NA NA NA NA
Experimental:T5-ABS1:T1 NA NA NA NA
ABS1:T6-ABS1:T1 NA NA NA NA
ABS2:T6-ABS1:T1 NA NA NA NA
Control:T6-ABS1:T1 NA NA NA NA
Experimental:T6-ABS1:T1 NA NA NA NA
Control:T1-ABS2:T1 NA NA NA NA
Experimental:T1-ABS2:T1 NA NA NA NA
ABS1:T4-ABS2:T1 NA NA NA NA
ABS2:T4-ABS2:T1 NA NA NA NA
Control:T4-ABS2:T1 NA NA NA NA
Experimental:T4-ABS2:T1 NA NA NA NA
ABS1:T5-ABS2:T1 NA NA NA NA
ABS2:T5-ABS2:T1 NA NA NA NA
Control:T5-ABS2:T1 NA NA NA NA
Experimental:T5-ABS2:T1 NA NA NA NA
ABS1:T6-ABS2:T1 NA NA NA NA
ABS2:T6-ABS2:T1 NA NA NA NA
Control:T6-ABS2:T1 NA NA NA NA
Experimental:T6-ABS2:T1 NA NA NA NA
Experimental:T1-Control:T1 -19.1 -61.526821 23.326821 0.9698630
ABS1:T4-Control:T1 -0.5 -48.805061 47.805061 1.0000000
ABS2:T4-Control:T1 -18.3 -60.726821 24.126821 0.9800159
Control:T4-Control:T1 -29.7 -72.126821 12.726821 0.4918565
Experimental:T4-Control:T1 NA NA NA NA
ABS1:T5-Control:T1 -52.9 -95.326821 -10.473179 0.0037057
ABS2:T5-Control:T1 -16.3 -58.726821 26.126821 0.9941168
Control:T5-Control:T1 -9.0 -53.721797 35.721797 0.9999994
Experimental:T5-Control:T1 NA NA NA NA
ABS1:T6-Control:T1 -53.3 -95.726821 -10.873179 0.0033305
ABS2:T6-Control:T1 -61.9 -104.326821 -19.473179 0.0003099
Control:T6-Control:T1 -17.7 -60.126821 24.726821 0.9857164
Experimental:T6-Control:T1 NA NA NA NA
ABS1:T4-Experimental:T1 18.6 -27.588473 64.788473 0.9901418
ABS2:T4-Experimental:T1 0.8 -39.200391 40.800391 1.0000000
Control:T4-Experimental:T1 -10.6 -50.600391 29.400391 0.9999546
Experimental:T4-Experimental:T1 NA NA NA NA
ABS1:T5-Experimental:T1 -33.8 -73.800391 6.200391 0.1922078
ABS2:T5-Experimental:T1 2.8 -37.200391 42.800391 1.0000000
Control:T5-Experimental:T1 10.1 -32.326821 52.526821 0.9999910
Experimental:T5-Experimental:T1 NA NA NA NA
ABS1:T6-Experimental:T1 -34.2 -74.200391 5.800391 0.1778735
ABS2:T6-Experimental:T1 -42.8 -82.800391 -2.799609 0.0250200
Control:T6-Experimental:T1 1.4 -38.600391 41.400391 1.0000000
Experimental:T6-Experimental:T1 NA NA NA NA
ABS2:T4-ABS1:T4 -17.8 -63.988473 28.388473 0.9939096
Control:T4-ABS1:T4 -29.2 -75.388473 16.988473 0.6667489
Experimental:T4-ABS1:T4 NA NA NA NA
ABS1:T5-ABS1:T4 -52.4 -98.588473 -6.211527 0.0127614
ABS2:T5-ABS1:T4 -15.8 -61.988473 30.388473 0.9985179
Control:T5-ABS1:T4 -8.5 -56.805061 39.805061 0.9999999
Experimental:T5-ABS1:T4 NA NA NA NA
ABS1:T6-ABS1:T4 -52.8 -98.988473 -6.611527 0.0116322
ABS2:T6-ABS1:T4 -61.4 -107.588473 -15.211527 0.0014377
Control:T6-ABS1:T4 -17.2 -63.388473 28.988473 0.9958800
Experimental:T6-ABS1:T4 NA NA NA NA
Control:T4-ABS2:T4 -11.4 -51.400391 28.600391 0.9998702
Experimental:T4-ABS2:T4 NA NA NA NA
ABS1:T5-ABS2:T4 -34.6 -74.600391 5.400391 0.1643710
ABS2:T5-ABS2:T4 2.0 -38.000391 42.000391 1.0000000
Control:T5-ABS2:T4 9.3 -33.126821 51.726821 0.9999975
Experimental:T5-ABS2:T4 NA NA NA NA
ABS1:T6-ABS2:T4 -35.0 -75.000391 5.000391 0.1516798
ABS2:T6-ABS2:T4 -43.6 -83.600391 -3.599609 0.0203729
Control:T6-ABS2:T4 0.6 -39.400391 40.600391 1.0000000
Experimental:T6-ABS2:T4 NA NA NA NA
Experimental:T4-Control:T4 NA NA NA NA
ABS1:T5-Control:T4 -23.2 -63.200391 16.800391 0.7901962
ABS2:T5-Control:T4 13.4 -26.600391 53.400391 0.9988609
Control:T5-Control:T4 20.7 -21.726821 63.126821 0.9387702
Experimental:T5-Control:T4 NA NA NA NA
ABS1:T6-Control:T4 -23.6 -63.600391 16.400391 0.7681937
ABS2:T6-Control:T4 -32.2 -72.200391 7.800391 0.2581287
Control:T6-Control:T4 12.0 -28.000391 52.000391 0.9997349
Experimental:T6-Control:T4 NA NA NA NA
ABS1:T5-Experimental:T4 NA NA NA NA
ABS2:T5-Experimental:T4 NA NA NA NA
Control:T5-Experimental:T4 NA NA NA NA
Experimental:T5-Experimental:T4 NA NA NA NA
ABS1:T6-Experimental:T4 NA NA NA NA
ABS2:T6-Experimental:T4 NA NA NA NA
Control:T6-Experimental:T4 NA NA NA NA
Experimental:T6-Experimental:T4 NA NA NA NA
ABS2:T5-ABS1:T5 36.6 -3.400391 76.600391 0.1085233
Control:T5-ABS1:T5 43.9 1.473179 86.326821 0.0356593
Experimental:T5-ABS1:T5 NA NA NA NA
ABS1:T6-ABS1:T5 -0.4 -40.400391 39.600391 1.0000000
ABS2:T6-ABS1:T5 -9.0 -49.000391 31.000391 0.9999963
Control:T6-ABS1:T5 35.2 -4.800391 75.200391 0.1456312
Experimental:T6-ABS1:T5 NA NA NA NA
Control:T5-ABS2:T5 7.3 -35.126821 49.726821 1.0000000
Experimental:T5-ABS2:T5 NA NA NA NA
ABS1:T6-ABS2:T5 -37.0 -77.000391 3.000391 0.0994925
ABS2:T6-ABS2:T5 -45.6 -85.600391 -5.599609 0.0120320
Control:T6-ABS2:T5 -1.4 -41.400391 38.600391 1.0000000
Experimental:T6-ABS2:T5 NA NA NA NA
Experimental:T5-Control:T5 NA NA NA NA
ABS1:T6-Control:T5 -44.3 -86.726821 -1.873179 0.0324703
ABS2:T6-Control:T5 -52.9 -95.326821 -10.473179 0.0037057
Control:T6-Control:T5 -8.7 -51.126821 33.726821 0.9999992
Experimental:T6-Control:T5 NA NA NA NA
ABS1:T6-Experimental:T5 NA NA NA NA
ABS2:T6-Experimental:T5 NA NA NA NA
Control:T6-Experimental:T5 NA NA NA NA
Experimental:T6-Experimental:T5 NA NA NA NA
ABS2:T6-ABS1:T6 -8.6 -48.600391 31.400391 0.9999982
Control:T6-ABS1:T6 35.6 -4.400391 75.600391 0.1341122
Experimental:T6-ABS1:T6 NA NA NA NA
Control:T6-ABS2:T6 44.2 4.199609 84.200391 0.0174282
Experimental:T6-ABS2:T6 NA NA NA NA
Experimental:T6-Control:T6 NA NA NA NA
head(p_df)
p_df %>% group_by(Time, Condition) %>% summarize(mean_asvs=mean(value), se_asv=sd(value)/sqrt(length(value)))
`summarise()` has grouped output by 'Time'. You can override using the `.groups` argument.
#replace NAs with 0s to get multcomp to run
tukey$`Condition:Time`[!complete.cases(tukey$`Condition:Time`),] <- 0
library(multcompView)
tukey.cld <-multcompLetters4(aov.trt, tukey)
print(tukey.cld)
$Condition
Control Experimental ABS2 ABS1
"a" "ab" "b" "b"
$Time
T0 T1 T4 T5 T6
"ab" "a" "a" "a" "b"
$`Condition:Time`
Control:T1 ABS1:T4 Control:T0 Control:T5 ABS2:T5 Control:T6 ABS2:T4 Experimental:T1
"a" "a" "a" "a" "ab" "ab" "ab" "ab"
Control:T4 ABS1:T5 ABS1:T6 ABS2:T6 ABS1:T0 ABS2:T0 tkl ABS1:T1
"abc" "bc" "bc" "c" "d" "e" "f" "g"
ABS2:T1 Experimental:T4 Experimental:T5 Experimental:T6
"h" "i" "j" "k"
tkl<-as.vector(tukey.cld$`Condition:Time`)
tkl.df<-data.frame(trt=tkl['Letters'])
tkl.df<-tkl.df %>% rownames_to_column() %>% separate(rowname,into=c('Condition','Time'))
head(tkl.df)
tkl.df$Condition<-factor(tkl.df$Condition, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))
tkl.df<-tkl.df[!tkl.df$Letters %in% c('d','e','f','g','h','i','j','k'),]
head(p_df)
p_df$Condition<-factor(p_df$Condition, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))
a<-ggplot(p_df, aes(y=value, x=Condition, fill=Condition)) +
geom_boxplot(alpha=0.5, outlier.shape = NA, color='black') +
geom_text(data=tkl.df, aes(y=10,x=Condition,label=Letters))+
scale_fill_manual(values=l2p16)+scale_color_manual(values=l2p16)+
geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+
theme_new + ylab('Nr. ASVs') + xlab('') + guides(fill=F, group=F,color=F) + ylim(c(0,140)) +
theme(axis.text.x = element_blank(),panel.grid = element_blank(),panel.grid.minor = element_blank(),rect = element_blank())
a
ASVs present in at least 90% of control samples at any relative abundance
# subset for control group
ps.control<-subset_samples(ps.red, Condition %in% c('Control'))
#removes taxa with 0 counts
ps.control <- prune_taxa(taxa_sums(ps.control) > 0, ps.control)
#transforms to rel. abundance
pseq.rel <- microbiome::transform(ps.control, "compositional")
plots across prevalance thresholds to help guide selection of the core
# Core with compositionals:
prevalences <- seq(.05, 1, .05)
detections <- round(10^seq(log10(1e-3), log10(.2), length = 10), 3)
# Also define gray color palette
gray <- gray(seq(0,1,length=5))
# Only include core taxa
pseq.core <- core(pseq.rel, detection =0, prevalence = 90/100)
p <- plot_core(pseq.core, plot.type = "heatmap",
colours = rev(brewer.pal(5, "Spectral")),
prevalences = prevalences,
detections = detections
) +
labs(x = "Detection Threshold")
Warning: The plot_core function is typically used with compositional
data. The data is not compositional. Make sure that you
intend to operate on non-compositional data.
print(p)
Subset entire relative abundance dataset with only core taxa
##identify core taxa
core.taxa <- core_members(pseq.rel, detection = 0, prevalence = 90/100)
ps.red.rel<-microbiome::transform(ps.red, "compositional")
ps.core<-subset_taxa(ps.red.rel, rownames(tax_table(ps.red.rel)) %in% core.taxa) #subsets relative abundance table from the ps.red sample dataset
ps.core
phyloseq-class experiment-level object
otu_table() OTU Table: [ 51 taxa and 54 samples ]
sample_data() Sample Data: [ 54 samples by 8 sample variables ]
tax_table() Taxonomy Table: [ 51 taxa by 7 taxonomic ranks ]
phy_tree() Phylogenetic Tree: [ 51 tips and 50 internal nodes ]
refseq() DNAStringSet: [ 51 reference sequences ]
51 Taxa make up the Aiptasia “core” as defined by prevalence in 90 % of all control samples across time points at any relative abundance
range(sample_sums(ps.core))
[1] 0.4119657 0.9966101
range(sample_sums(subset_samples(ps.core, Condition=='Experimental')))
[1] 0.4119657 0.7929611
range(sample_sums(subset_samples(ps.core, Condition=='ABS1')))
[1] 0.7805256 0.9661604
range(sample_sums(subset_samples(ps.core, Condition=='ABS2')))
[1] 0.7163705 0.9966101
microbiome::plot_composition(ps.core,
plot.type = "barplot", sample.sort = "neatmap",group_by = 'Condition')
core.taxa.df<-data.frame(tax_table(ps.core)) %>% arrange(Phylum, Class)
core.taxa.df
core.taxa.df %>% group_by(Class) %>% summarize(nr=n())
# ordinate
set.seed(99)
ord<-ordinate(ps.core,"PCoA",distance = "bray")
ord.df<-plot_ordination(ps.core,ord,justDF=T)
fb<-ggplot(ord.df, aes(x=Axis.1, y =Axis.2, color=Condition, shape = Time, Fill=Condition)) +
geom_point(size=4) +
scale_shape_manual(values=shpsc) +
scale_color_manual(values=l2p16) + scale_fill_manual(values=l2p16) +
theme_new #+ guides(shape=F, color=F)
fb
betadispersion
#run adonis
iDist <- phyloseq::distance(ps.core, method = "bray") # note, must call phyloseq specifically if deseq2 is also loaded
pn_df <- data.frame(sample_data(ps.core))
#betadispersion
samps<-data.frame(sample_data(ps.red))
#betdisper
bd<-betadisper(iDist, paste(samps$Time,samps$Condition))
bd.df<-data.frame(distances=bd$distances, group1=bd$group,group2=bd$group)
bd.df<-bd.df %>% separate(group1,into = c('Time','Treatment'))
bd.df$Treatment<-factor(bd.df$Treatment, ordered=T, levels=c('Control','Experimental','ABS1','ABS2'))
c<-ggplot(bd.df, aes(group=group2, y=distances, x=Treatment, fill=Treatment)) +
geom_boxplot(alpha=0.5, outlier.shape = NA,color='black') +
scale_fill_manual(values=l2p16)+ scale_color_manual(values=l2p16)+
geom_jitter(width=0.1,height=0,size=1, alpha=0.5) + facet_grid(~Time,scales='free',space='free')+
theme_new + ylab('Distance to centriod') + xlab('') + guides(fill=F, group=F, color=F) +
theme(axis.text.x = element_blank(),panel.grid = element_blank(),plot.background = element_blank(),panel.grid.minor = element_blank(),rect = element_blank())
c
stats
# Adonis test
mod<-adonis2(iDist ~ Condition*Time, data = pn_df)
mod
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = iDist ~ Condition * Time, data = pn_df)
Df SumOfSqs R2 F Pr(>F)
Condition 3 2.5726 0.18752 7.7992 0.001 ***
Time 4 4.5488 0.33156 10.3425 0.001 ***
Condition:Time 4 1.9799 0.14431 4.5016 0.001 ***
Residual 42 4.6181 0.33661
Total 53 13.7194 1.00000
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
plot(mod)
#betadispersion
aov<-anova(bd)
aov
Analysis of Variance Table
Response: Distances
Df Sum Sq Mean Sq F value Pr(>F)
Groups 11 0.44793 0.040721 1.2764 0.2711
Residuals 42 1.33996 0.031904
bd
Homogeneity of multivariate dispersions
Call: betadisper(d = iDist, group = paste(samps$Time, samps$Condition))
No. of Positive Eigenvalues: 36
No. of Negative Eigenvalues: 17
Average distance to median:
T0 Control T1 Control T1 Experimental T4 ABS1 T4 ABS2 T4 Control T5 ABS1 T5 ABS2
0.1688 0.1995 0.2540 0.1112 0.1887 0.2997 0.4210 0.3990
T5 Control T6 ABS1 T6 ABS2 T6 Control
0.2133 0.1486 0.2152 0.2928
Eigenvalues for PCoA axes:
(Showing 8 of 53 eigenvalues)
PCoA1 PCoA2 PCoA3 PCoA4 PCoA5 PCoA6 PCoA7 PCoA8
4.6413 2.5189 1.6697 0.9451 0.7118 0.5287 0.4723 0.3704
tukey.bd<-TukeyHSD(bd, ordered=T)
tukey.bd
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = distances ~ group, data = df)
$group
diff lwr upr p adj
T6 ABS1-T4 ABS1 0.037406752 -0.4136519 0.4884654 1.0000000
T0 Control-T4 ABS1 0.057603553 -0.4466953 0.5619025 0.9999997
T4 ABS2-T4 ABS1 0.077534137 -0.3735245 0.5285928 0.9999773
T1 Control-T4 ABS1 0.088337882 -0.3833905 0.5600663 0.9999466
T5 Control-T4 ABS1 0.102137167 -0.3695913 0.5738656 0.9997794
T6 ABS2-T4 ABS1 0.104036920 -0.3470217 0.5550956 0.9995975
T1 Experimental-T4 ABS1 0.142810858 -0.3082478 0.5938695 0.9933197
T6 Control-T4 ABS1 0.181585017 -0.2694736 0.6326437 0.9587465
T4 Control-T4 ABS1 0.188454688 -0.2626040 0.6395133 0.9469959
T5 ABS2-T4 ABS1 0.287766418 -0.1632922 0.7388251 0.5545898
T5 ABS1-T4 ABS1 0.309798021 -0.1412606 0.7608567 0.4439119
T0 Control-T6 ABS1 0.020196802 -0.4308618 0.4712554 1.0000000
T4 ABS2-T6 ABS1 0.040127385 -0.3505009 0.4307556 0.9999999
T1 Control-T6 ABS1 0.050931130 -0.3633927 0.4652550 0.9999993
T5 Control-T6 ABS1 0.064730416 -0.3495934 0.4790542 0.9999914
T6 ABS2-T6 ABS1 0.066630169 -0.3239981 0.4572584 0.9999790
T1 Experimental-T6 ABS1 0.105404106 -0.2852241 0.4960324 0.9983034
T6 Control-T6 ABS1 0.144178266 -0.2464500 0.5348065 0.9778134
T4 Control-T6 ABS1 0.151047936 -0.2395803 0.5416762 0.9688830
T5 ABS2-T6 ABS1 0.250359666 -0.1402686 0.6409879 0.5478056
T5 ABS1-T6 ABS1 0.272391269 -0.1182370 0.6630195 0.4211407
T4 ABS2-T0 Control 0.019930583 -0.4311281 0.4709892 1.0000000
T1 Control-T0 Control 0.030734329 -0.4409941 0.5024628 1.0000000
T5 Control-T0 Control 0.044533614 -0.4271948 0.5162620 1.0000000
T6 ABS2-T0 Control 0.046433367 -0.4046253 0.4974920 0.9999999
T1 Experimental-T0 Control 0.085207305 -0.3658513 0.5362660 0.9999417
T6 Control-T0 Control 0.123981464 -0.3270772 0.5750401 0.9980013
T4 Control-T0 Control 0.130851135 -0.3202075 0.5819098 0.9968010
T5 ABS2-T0 Control 0.230162864 -0.2208958 0.6812215 0.8272338
T5 ABS1-T0 Control 0.252194468 -0.1988642 0.7032531 0.7324375
T1 Control-T4 ABS2 0.010803745 -0.4035201 0.4251276 1.0000000
T5 Control-T4 ABS2 0.024603031 -0.3897208 0.4389269 1.0000000
T6 ABS2-T4 ABS2 0.026502784 -0.3641255 0.4171310 1.0000000
T1 Experimental-T4 ABS2 0.065276721 -0.3253515 0.4559050 0.9999830
T6 Control-T4 ABS2 0.104050881 -0.2865774 0.4946791 0.9984881
T4 Control-T4 ABS2 0.110920551 -0.2797077 0.5015488 0.9973420
T5 ABS2-T4 ABS2 0.210232281 -0.1803960 0.6008605 0.7752301
T5 ABS1-T4 ABS2 0.232263884 -0.1583644 0.6228921 0.6544875
T5 Control-T1 Control 0.013799285 -0.4229364 0.4505349 1.0000000
T6 ABS2-T1 Control 0.015699038 -0.3986248 0.4300229 1.0000000
T1 Experimental-T1 Control 0.054472976 -0.3598508 0.4687968 0.9999986
T6 Control-T1 Control 0.093247135 -0.3210767 0.5075710 0.9996809
T4 Control-T1 Control 0.100116806 -0.3142070 0.5144406 0.9993781
T5 ABS2-T1 Control 0.199428536 -0.2148953 0.6137524 0.8738175
T5 ABS1-T1 Control 0.221460139 -0.1928637 0.6357840 0.7824422
T6 ABS2-T5 Control 0.001899753 -0.4124241 0.4162236 1.0000000
T1 Experimental-T5 Control 0.040673691 -0.3736501 0.4549975 0.9999999
T6 Control-T5 Control 0.079447850 -0.3348760 0.4937717 0.9999324
T4 Control-T5 Control 0.086317521 -0.3280063 0.5006413 0.9998478
T5 ABS2-T5 Control 0.185629251 -0.2286946 0.5999531 0.9172050
T5 ABS1-T5 Control 0.207660854 -0.2066630 0.6219847 0.8426829
T1 Experimental-T6 ABS2 0.038773938 -0.3518543 0.4294022 0.9999999
T6 Control-T6 ABS2 0.077548097 -0.3130802 0.4681763 0.9999049
T4 Control-T6 ABS2 0.084417768 -0.3062105 0.4750460 0.9997834
T5 ABS2-T6 ABS2 0.183729497 -0.2068988 0.5743577 0.8893728
T5 ABS1-T6 ABS2 0.205761101 -0.1848671 0.5963893 0.7973195
T6 Control-T1 Experimental 0.038774159 -0.3518541 0.4294024 0.9999999
T4 Control-T1 Experimental 0.045643830 -0.3449844 0.4362721 0.9999996
T5 ABS2-T1 Experimental 0.144955560 -0.2456727 0.5355838 0.9769141
T5 ABS1-T1 Experimental 0.166987163 -0.2236411 0.5576154 0.9384202
T4 Control-T6 Control 0.006869671 -0.3837586 0.3974979 1.0000000
T5 ABS2-T6 Control 0.106181400 -0.2844468 0.4968096 0.9981889
T5 ABS1-T6 Control 0.128213004 -0.2624152 0.5188413 0.9910467
T5 ABS2-T4 Control 0.099311730 -0.2913165 0.4899400 0.9990076
T5 ABS1-T4 Control 0.121343333 -0.2692849 0.5119716 0.9942930
T5 ABS1-T5 ABS2 0.022031603 -0.3685966 0.4126599 1.0000000
#replace NAs with 0s to get multcomp to run
tukey$`Condition:Time`[!complete.cases(tukey$`Condition:Time`),] <- 0
library(multcompView)
tukey.cld <-multcompLetters4(aov,comp = tukey.bd)
print(tukey.cld)
multcompLetters(tukey.bd$group[,4])
#Fig 2D plot distance tree
tr<-plot_tree(ps.core, 'treeonly',label.tips="taxa_names",nodelabf=nodeplotblank,ladderize="left", justify = "yes please")
Warning: Some columns are a multi-column type (such as a matrix column): [1, 2, 3, 4, 5, 6, 7]. setDT will retain these columns as-is but subsequent operations like grouping and joining may fail. Please consider as.data.table() instead which will create a new column for each embedded column.Warning: Some columns are a multi-column type (such as a matrix column): [1, 2, 3, 4, 5, 6, 7]. setDT will retain these columns as-is but subsequent operations like grouping and joining may fail. Please consider as.data.table() instead which will create a new column for each embedded column.
asv.order<-data.frame(tr$data) %>% arrange(y)
asv.order<-na.omit(as.vector(asv.order$OTU))
asv.order
[1] "ASV1" "ASV5" "ASV56" "ASV3" "ASV19" "ASV6" "ASV86" "ASV4" "ASV9" "ASV25" "ASV24" "ASV21" "ASV61" "ASV110"
[15] "ASV44" "ASV22" "ASV35" "ASV39" "ASV16" "ASV37" "ASV67" "ASV13" "ASV152" "ASV51" "ASV95" "ASV15" "ASV81" "ASV27"
[29] "ASV47" "ASV34" "ASV8" "ASV90" "ASV75" "ASV26" "ASV12" "ASV14" "ASV32" "ASV79" "ASV87" "ASV2" "ASV49" "ASV91"
[43] "ASV40" "ASV101" "ASV43" "ASV38" "ASV59" "ASV7" "ASV20" "ASV17" "ASV33"
attr(,"na.action")
[1] 4 6 8 10 11 13 17 19 20 21 23 25 27 31 32 34 36 37 39 40 44 46 48 50 52 54 55 56 58 60 63 65 67 70 71 73 76 77 79 81 84 85
[43] 87 89 91 93 95 97 99
attr(,"class")
[1] "omit"
d<-ggtree::ggtree(ps.core) + ggtree::geom_tiplab(align=T, hjust=-.3) + ggtree::hexpand(.01)
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
d
Error in `geom_segment2()`:
! Problem while converting geom to grob.
ℹ Error occurred in the 3rd layer.
Caused by error in `check.length()`:
! 'gpar' element 'lwd' must not be length 0
Backtrace:
1. base (local) `<fn>`(x)
2. ggplot2:::print.ggplot(x)
4. ggplot2:::ggplot_gtable.ggplot_built(data)
5. ggplot2:::by_layer(...)
12. ggplot2 (local) f(l = layers[[i]], d = data[[i]])
...
24. grid::gpar(...)
25. grid:::validGP(list(...))
26. grid (local) numnotnull("lwd")
27. grid (local) check.length(gparname)
28. base::stop(...)
library(ggtree)
ggtree v3.4.4 For help: https://yulab-smu.top/treedata-book/
If you use the ggtree package suite in published research, please cite the appropriate paper(s):
Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and
annotation of phylogenetic trees with their covariates and other associated data. Methods in Ecology and Evolution.
2017, 8(1):28-36. doi:10.1111/2041-210X.12628
LG Wang, TTY Lam, S Xu, Z Dai, L Zhou, T Feng, P Guo, CW Dunn, BR Jones, T Bradley, H Zhu, Y Guan, Y Jiang, G Yu.
treeio: an R package for phylogenetic tree input and output with richly annotated and associated data. Molecular
Biology and Evolution. 2020, 37(2):599-603. doi: 10.1093/molbev/msz240
Guangchuang Yu, Tommy Tsan-Yuk Lam, Huachen Zhu, Yi Guan. Two methods for mapping and visualizing associated data on
phylogeny using ggtree. Molecular Biology and Evolution. 2018, 35(12):3041-3043. doi:10.1093/molbev/msy194
Attaching package: ‘ggtree’
The following object is masked from ‘package:Biostrings’:
collapse
The following object is masked from ‘package:IRanges’:
collapse
The following object is masked from ‘package:S4Vectors’:
expand
The following object is masked from ‘package:ape’:
rotate
The following object is masked from ‘package:tidyr’:
expand
ggtree(ps.core) + geom_tiplab(align=F, hjust=-.3, aes(label=Class))
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
# Bubble plot of relative abundances in ps.core
# Include only 1 control sample per time point to reduce complexity of plot
ps.core.control<-subset_samples(ps.core,Condition=='Control') %>% subset_samples(Sample_names %in%c('SAM_133_S1','SAM_157_S10','SAM_347_S56','SAM_364_S73','SAM_409_S89'))
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
ps.core.other<-subset_samples(ps.core,!Condition=='Control')
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
#merge
ps.core.new<-merge_phyloseq(ps.core.control,ps.core.other)
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
# place samples into correct order
control.samps<-data.frame(sample_data(ps.core.control))$Sample_names
other.samps<-data.frame(sample_data(ps.core.other))$Sample_names
samps.order<-c(control.samps,other.samps)
#Add Grouping variable for plot
samps.new<-data.frame(sample_data(ps.core.new)) %>% mutate(PlotGroup=case_when(Condition=='Control' ~ 'ASW',
Condition=='Experimental' ~ 'FASW 33 d',
Condition=='ABS1' & Time=='T4' ~ "ABS1 55 d",
Condition=='ABS1' & Time=='T5' ~ "ABS1 61 d",
Condition=='ABS1' & Time=='T6' ~ "ABS1 76 d",
Condition=='ABS2' & Time=='T4' ~ "ABS2 55 d",
Condition=='ABS2' & Time=='T5' ~ "ABS2 61 d",
Condition=='ABS2' & Time=='T6' ~ "ABS2 76 d"))
sample_data(ps.core.new)<-samps.new
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
Also defined by ‘tidytree’
#make color palette for core taxa via family level classification
df.taxa<-data.frame(ASV=rownames(tax_table(ps.core.new)),tax_table(ps.core.new)) %>% arrange(Phylum, Class, Order, Family, Genus)
df.taxa<-data.frame(ASV=rownames(tax_table(ps.core.new)),tax_table(ps.core.new)) %>% arrange(factor(ASV, levels=asv.order))
df.taxa$Family2<-ifelse(is.na(df.taxa$Family),paste(df.taxa$Class,'-unknown',sep=''),df.taxa$Family) #fix na family level classificaiton
df.taxa$Family2<-ifelse(df.taxa$Family2=='NA-unknown',paste(df.taxa$Phylum,'-unknown',sep=''),df.taxa$Family2) #fix na family level classificaiton
df.taxa
#set color pallet based on taxa order in tree, reds gammas, yellows, alphas, purples
clr<-c('#FCD9DF','#38050D','#F7A1AF','#F05670','#EB1E40','#CE1231','#960D24', #"Alteromonadaceae","Gammaproteobacteria-unknown","Vibrionaceae","Halomonadaceae","Pseudomonadaceae","Oleiphilaceae",Spongiibacteraceae
'#D46021',# Akkermansiaceae
'mediumpurple4','mediumpurple1', #"Phycisphaerae-unknown","Pirellulaceae"
'#9BC0AF', #Nannocystaceae
'#187B4D', #Mycobacteriaceae
'#093149','#49AEE9','#A4D7F4',#"Bacteroidia-unknown","Saprospiraceae","Flavobacteriaceae"
'#F9F5DC',"#F0E7A8",'#F7EB8D',"#E7DA73",#"Emcibacteraceae","Hyphomonadaceae","GCA-2696645","Sphingomonadaceae"
"gray",#"Proteobacteria-unknown"
'#F3E04B','#574F0F','#F0DA32','#E4CB11','#BEA90E','#85760A'#"Devosiaceae","Stappiaceae","Alphaproteobacteria-unknown", "Cohaesibacteraceae"
#"Rhodobacteraceae","Rhizobiaceae"
)
names(clr)<-unique(df.taxa$Family2)
make bubble plot
# get relavent dataframes
md<-data.frame(sample_data(ps.core.new))
df.asv<-data.frame(Sample_names=rownames(otu_table(ps.core.new)),otu_table(ps.core.new))
#make long dataframe for plotting
df.asv.long<-pivot_longer(df.asv,cols=!Sample_names,names_to='ASV',values_to = 'RelAbund') #melt dataframe
df.relabund<-df.asv.long %>% left_join(y=md, by='Sample_names') #add in metadata
df.relabund<-df.relabund %>% left_join(y=df.taxa, by='ASV') #add in taxa info
df.relabund<-df.relabund[df.relabund$RelAbund>0,]
#re-order ASVs
df.relabund$ASV<-factor(df.relabund$ASV, ordered=T, levels=asv.order)
#reoder samples
df.relabund$Sample_names<-factor(df.relabund$Sample_names, ordered=T, levels=samps.order)
#re-order treatments
unique(df.relabund$PlotGroup)
[1] "ASW" "FASW 33 d" "ABS1 55 d" "ABS2 55 d" "ABS1 61 d" "ABS2 61 d" "ABS1 76 d" "ABS2 76 d"
df.relabund$PlotGroup<-factor(df.relabund$PlotGroup, ordered=T, levels=c("ASW","FASW 33 d", "ABS1 55 d", "ABS2 55 d", "ABS1 61 d", "ABS2 61 d", "ABS1 76 d","ABS2 76 d"))
head(df.relabund)
e<-ggplot(df.relabund, aes(x = Sample_names, y = ASV, fill = Family2)) +
geom_point(aes(size = RelAbund*10), alpha = 0.8, shape = 21) +
scale_size_continuous(limits = c(0, 70), range = c(0,20), breaks = c(0.01, 5, 10, 20, 40, 80)) +
theme_new + xlab('') +ylab('') +
facet_grid(~PlotGroup, scales='free',space='free')+
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
labs(size = "Relative Abundance", color = "Family") +
scale_fill_manual(values=clr) + guides(fill=F)
e
#ggsave('maggie_plots/bubbles-with-legend.pdf')
put all together
library(patchwork)
layout3<-"AABBBB
CCBBBB
DEEEEE
DEEEEE
DEEEEE
"
e2<-e+theme(panel.grid = element_blank(),rect = element_blank(),panel.grid.minor = element_blank(),panel.grid.major = element_blank())
e2
fb2<-fb + guides(color="none", fill="none", shape="none")
a + fb2 + c + d + e2 + plot_layout(design=layout3)
ggsave('Figures/Figure2-new.pdf', width=11,height=10)
#save core taxa as .csv
df.taxa<-df.taxa %>% unite('string',c(Phylum,Class,Order,Family,Genus,Species), sep='_;',na.rm=F, remove=F)
head(df.taxa)
write.csv(df.taxa, 'Figures/core.taxa.csv')
Figure S1
samps<-data.frame(sample_data(ps.core))
samps$Condition<- factor(samps$Condition, ordered=T, levels=c('Control', 'Experimental','ABS1','ABS2'))
samps<-samps %>% arrange(Condition,Time)
#make new grouping factor
samps$grp<-paste(samps$Condition, samps$Time)
samps$grp<-factor(samps$grp, ordered = T, levels=c("Control T0","Control T1","Control T4", "Control T5","Control T6","Experimental T1", "ABS1 T4","ABS1 T5","ABS1 T6","ABS2 T4", "ABS2 T5","ABS2 T6"))
sample_data(ps.core)<-samps
ps.core %>%
aggregate_taxa(level='Family') %>%
plot_composition(plot.type='barplot',group_by = 'grp',sample.sort = 'neatmap') +
scale_fill_manual(values=clr)+
scale_y_continuous(label = scales::percent) +
theme_new + guides(fill=F)+
# Removes sample names and ticks
theme( axis.text.x = element_blank(), axis.ticks.x=element_blank(), plot.background = element_blank())
ggsave('Figures/figS1.pdf',height=6,width=12)